##########################################################################
####################        ## ANALYSES FULL ##       ####################
#####################           Max Joosten         ######################
#####################   Last updated: 14.06.2023    ######################
##########################################################################

# Read packages
library(stargazer)
library(vtable)
library(tidyverse)
library(sandwich)
library(sjPlot)
library(ggeffects)

# set working directory
setwd("your path here")
load("df_complete_v2_check.RData")

# Some minor changes
df$yeardistance <- df$year - df$orig.year
df$year <- as.factor(df$year)
df$country <- as.factor(df$country)
df <- df %>% mutate(spend_polarizing = (spend_health + spend_pension + spend_unemp)/3)

############################################################################
# Main text Table 1: Average spending preferences across income groups
df %>% filter(!is.na(cut_spend) & !is.na(l_ip_per414_mean) & !is.na(l_ip_per414_consensus) & !is.na(l_ip_per414) & !is.na(country) & !is.na(year) & !is.na(unemp) & !is.na(growth) & !is.na(age) & !is.na(gender)) %>%
  group_by(inc.group) %>%
  summarise(mean_cut = mean(cut_spend, na.rm = T), 
            mean_spend = mean(spend_polarizing, na.rm = T)) %>%
  filter(!is.na(inc.group))

############################################################################
# Main text Table 2
m1.0 <- lm(cut_spend ~   l_ip_per414_mean * l_ip_per414_consensus + l_ip_per414 + country + year + loggdp + unemp + growth + age + gender, weights = weight, data = df)
m1 <-   lm(cut_spend ~   l_ip_per414_mean * l_ip_per414_consensus + l_ip_per414 + country + year + loggdp + unemp + growth + age + gender, weights = weight, data = subset(df, inc.group==0))
m2 <-   lm(cut_spend ~   l_ip_per414_mean * l_ip_per414_consensus + l_ip_per414 + country + year + loggdp + unemp + growth + age + gender, weights = weight, data = subset(df, inc.group==1))
m3 <-   lm(cut_spend ~   l_ip_per414_mean * l_ip_per414_consensus + l_ip_per414 + country + year + loggdp + unemp + growth + age + gender, weights = weight, data = subset(df, inc.group==2))

## Robust SE: Baseline model
r1.0 <-  sqrt(diag(vcovHC(m1.0, type = "HC1")))
r1 <-  sqrt(diag(vcovHC(m1, type = "HC1")))
r2 <-  sqrt(diag(vcovHC(m2, type = "HC1")))
r3 <-  sqrt(diag(vcovHC(m3, type = "HC1")))

# Regression
stargazer(m1.0, m1, m2, m3,
          se = list(r1.0, r1, r2, r3),
          keep.stat=c("n", "adj.rsq"),
          type = "text", 
          # out = "../05_figures/v3/main1.html",
          omit = c("country", "year"),
          covariate.labels = c("Economic orthodoxy (average, t-1)", "Economic orthodoxy (consensus, t-1)", "Economic orthodoxy (party, t-1)", "GDP (logged) (t)", "Unemployment rate (t)", "Growth (t)", "Age", "Female", "Economic orthodoxy average * consensus"),
          dep.var.labels = c("Preference for spending cuts"),
          add.lines=list(c("Income group", "All",  "L", "M", "H"),
                         c("Country-years", "60", "60", "60", "60"),
                         c("Fixed effects", "C&Y", "C&Y", "C&Y", "C&Y", "C&Y", "C&Y", "C&Y", "C&Y", "C&Y")))


# Remove all models, only keep the dataframe
rm(list=setdiff(ls(), "df"))

############################################################################
# Main text Table 3
m1.0 <- lm(spend_polarizing ~   l_ip_per414_mean * l_ip_per414_consensus + l_ip_per414 + country + year + loggdp + unemp + growth + age + gender, weights = weight, data = df)
m1 <-   lm(spend_polarizing ~   l_ip_per414_mean * l_ip_per414_consensus + l_ip_per414 + country + year + loggdp + unemp + growth + age + gender, weights = weight, data = subset(df, inc.group==0))
m2 <-   lm(spend_polarizing ~   l_ip_per414_mean * l_ip_per414_consensus + l_ip_per414 + country + year + loggdp + unemp + growth + age + gender, weights = weight, data = subset(df, inc.group==1))
m3 <-   lm(spend_polarizing ~   l_ip_per414_mean * l_ip_per414_consensus + l_ip_per414 + country + year + loggdp + unemp + growth + age + gender, weights = weight, data = subset(df, inc.group==2))

## Robust SE: Baseline model
r1.0 <-  sqrt(diag(vcovHC(m1.0, type = "HC1")))
r1 <-  sqrt(diag(vcovHC(m1, type = "HC1")))
r2 <-  sqrt(diag(vcovHC(m2, type = "HC1")))
r3 <-  sqrt(diag(vcovHC(m3, type = "HC1")))

# Regression
stargazer(m1.0, m1, m2, m3,
          se = list(r1.0, r1, r2, r3),
          keep.stat=c("n", "adj.rsq"),
          type = "text", 
          #out = "../05_figures/v3/main2.html",
          omit = c("country", "year"),
          covariate.labels = c("Economic orthodoxy (average, t-1)", "Economic orthodoxy (consensus, t-1)", "Economic orthodoxy (party, t-1)", "GDP (logged) (t)", "Unemployment rate (t)", "Growth (t)", "Age", "Female", "Economic orthodoxy average * consensus"),
          dep.var.labels = c("Preference for spending in pensions, health and unemployment"),
          add.lines=list(c("Income group", "All",  "L", "M", "H"),
                         c("Country-years", "60", "60", "60", "60"),
                         c("Fixed effects", "C&Y", "C&Y", "C&Y", "C&Y", "C&Y", "C&Y", "C&Y", "C&Y", "C&Y")))

############################################################################################################
################################# Appendix A: A.	Descriptive information ##################################
############################################################################################################

# Appendix A1: summary statistics
labs <- data.frame(name_df = c("country", "year", "party", "inc.group", "gender", "age",
                               "cut_spend",
                               "spend_env", "spend_health", "spend_law", "spend_edu", "spend_def", "spend_pension", "spend_unemp", "spend_culture", "spend_total", 
                               "tax_L", "tax_M", "tax_H", "tax_total",
                               "l_ip_per414_mean", "l_ip_per414_consensus",
                               "l_ip_per409_mean", "l_ip_per409_consensus",
                               "growth", "unemp", "loggdp"),

label_df = c("Country", "Year", "Party", "Income group", "Gender", "Age",
             "Pref. Spending cuts",
             "Pref. spend env.", "Pref. spend health", "Pref. spend law", "Pref. spend edu.", "Pref. spend def.", "Pref. spend pension", "Pref. spend unemp.", "Pref. spend cult.", "Pref. spend total", 
             "Pref. tax (Low)", "Pref. tax (M)", "Pref. tax (H)", "Pref. tax (total)",
             "Economic orthodoxy average (t-1, interpolated, standardized)", "Economic orthodoxy consensus (t-1, interpolated, standardized)",
             "Keynesian demand management (t-1, interpolated, standardized)", "Keynesian demand management (t-1, interpolated, standardized)",
             "GDP growth", "Unemployment rate", "GDP (logged)"
             
))

dft <- df %>% filter(!is.na(inc.group) & !is.na(party))

st(dft, vars = c("year",
                 "inc.group", "gender", "age",
                 "cut_spend",
                 "spend_total", "spend_env", "spend_health", "spend_law", "spend_edu", "spend_def", "spend_pension", "spend_unemp", "spend_culture",
                 "tax_total", "tax_L", "tax_M", "tax_H",
                 "l_ip_per414_mean", "l_ip_per414_consensus",
                 "l_ip_per409_mean", "l_ip_per409_consensus",
                 "growth", "unemp", "loggdp"),
   summ = list(c('mean(x)','sd(x)','min(x)','max(x)', 'notNA(x)')),
   summ.names = list(c("Mean", "SD", "Min.", "Max.", "N")),
   labels = labs,
   file = "../05_figures/v3/appendix_A1.html")

# Appendix A2: Country-year average party attention to fiscal strategy, consensus and grouped preferences by income group
dtable <- dft %>% group_by(country, year, l_ip_per414_mean, l_ip_per414_consensus, l_ip_per409_mean, l_ip_per409_consensus, inc.group) %>% summarise(mean = mean(cut_spend, na.rm = T)) %>%
  pivot_wider(names_from = c("inc.group"), values_from = c("mean")) %>% 
  rename("Low" = `0`,
         "Middle" = `1`,
         "High" = `2`,
         "Country" = country,
         "Year" = year,
         "Economic orthodoxy average" = l_ip_per414_mean, 
         "Economic orthodoxy consensus" = l_ip_per414_consensus,
         "Keynesian demand management average" = l_ip_per409_mean, 
         "Keynesian demand management consensus" = l_ip_per409_consensus) # %>% select(!`NA`)

dtable$`Economic orthodoxy average` <- round(dtable$`Economic orthodoxy average`, digits = 3)
dtable$`Economic orthodoxy consensus` <- round(dtable$`Economic orthodoxy consensus`, digits = 3)
dtable$`Keynesian demand management average` <- round(dtable$`Keynesian demand management average`, digits = 3)
dtable$`Keynesian demand management consensus` <- round(dtable$`Keynesian demand management consensus`, digits = 3)
dtable$Low <- round(dtable$Low, digits = 2)
dtable$Middle <- round(dtable$Middle, digits = 2)
dtable$High <- round(dtable$High, digits = 2)

write.table(dtable, file = "../05_figures/v3/appendix_A2.txt", sep = ";", quote = FALSE, row.names = F)

############################################################################
# Appendix A3: Number of respondents per income group per country-year
dtable2 <- dft %>% 
  group_by(country, year, inc.group) %>% 
  summarize(n = n()) %>% 
  pivot_wider(names_from = c("inc.group"), values_from = c("n"))

write.table(dtable2, file = "../05_figures/v3/appendix_A3.txt", sep = ";", quote = FALSE, row.names = F)

############################################################################################################
###################### Appendix B: Results by individual spending and taxation items #######################
############################################################################################################

# Table B1: Differences across spending areas:
df %>% 
  group_by(inc.group) %>%
  summarise(mean_spend_env = mean(spend_env, na.rm = T),
            mean_spend_health = mean(spend_health, na.rm = T), 
            mean_spend_law = mean(spend_law, na.rm = T),
            mean_spend_edu = mean(spend_edu, na.rm = T),
            mean_spend_def = mean(spend_def, na.rm = T),
            mean_spend_pension = mean(spend_pension, na.rm = T),
            mean_spend_unemp = mean(spend_unemp, na.rm = T),
            mean_spend_culture = mean(spend_culture, na.rm = T),
            
            mean_tax_L = mean(tax_L, na.rm = T),
            mean_tax_M = mean(tax_M, na.rm = T),
            mean_tax_H = mean(tax_H, na.rm = T),
  ) %>%
  filter(!is.na(inc.group))

# Appendix B2: spend pref low
m1 <-   lm(spend_env ~   l_ip_per414_mean * l_ip_per414_consensus + l_ip_per414 + as.factor(country) + as.factor(year) + loggdp + unemp + growth + age + gender, weights = weight, data = subset(df, inc.group==0))
m2 <-   lm(spend_health ~   l_ip_per414_mean * l_ip_per414_consensus + l_ip_per414 + as.factor(country) + as.factor(year) + loggdp + unemp + growth + age + gender, weights = weight, data = subset(df, inc.group==0))
m3 <-   lm(spend_law ~   l_ip_per414_mean * l_ip_per414_consensus + l_ip_per414 + as.factor(country) + as.factor(year) + loggdp + unemp + growth + age + gender, weights = weight, data = subset(df, inc.group==0))
m4 <-   lm(spend_edu ~   l_ip_per414_mean * l_ip_per414_consensus + l_ip_per414 + as.factor(country) + as.factor(year) + loggdp + unemp + growth + age + gender, weights = weight, data = subset(df, inc.group==0))
m5 <-   lm(spend_def ~   l_ip_per414_mean * l_ip_per414_consensus + l_ip_per414 + as.factor(country) + as.factor(year) + loggdp + unemp + growth + age + gender, weights = weight, data = subset(df, inc.group==0))
m6 <-   lm(spend_pension ~   l_ip_per414_mean * l_ip_per414_consensus + l_ip_per414 + as.factor(country) + as.factor(year) + loggdp + unemp + growth + age + gender, weights = weight, data = subset(df, inc.group==0))
m7 <-   lm(spend_unemp ~   l_ip_per414_mean * l_ip_per414_consensus + l_ip_per414 + as.factor(country) + as.factor(year) + loggdp + unemp + growth + age + gender, weights = weight, data = subset(df, inc.group==0))
m8 <-   lm(spend_culture ~   l_ip_per414_mean * l_ip_per414_consensus + l_ip_per414 + as.factor(country) + as.factor(year) + loggdp + unemp + growth + age + gender, weights = weight, data = subset(df, inc.group==0))
m9 <-   lm(spend_total ~   l_ip_per414_mean * l_ip_per414_consensus + l_ip_per414 + as.factor(country) + as.factor(year) + loggdp + unemp + growth + age + gender, weights = weight, data = subset(df, inc.group==0))

## Robust SE: Baseline model
r1 <-  sqrt(diag(vcovHC(m1, type = "HC1")))
r2 <-  sqrt(diag(vcovHC(m2, type = "HC1")))
r3 <-  sqrt(diag(vcovHC(m3, type = "HC1")))
r4 <-  sqrt(diag(vcovHC(m4, type = "HC1")))
r5 <-  sqrt(diag(vcovHC(m5, type = "HC1")))
r6 <-  sqrt(diag(vcovHC(m6, type = "HC1")))
r7 <-  sqrt(diag(vcovHC(m7, type = "HC1")))
r8 <-  sqrt(diag(vcovHC(m8, type = "HC1")))
r9 <-  sqrt(diag(vcovHC(m9, type = "HC1")))
# Regression
stargazer(m1, m2, m3, m4, m5, m6, m7, m8, m9,
          se = list(r1, r2, r3, r4, r5, r6, r7, r8, r9),
          keep.stat=c("n", "adj.rsq"),
          type = "text", 
          #out = "../05_figures/v3/appendix.ortho.lowspend.html",
          title = "Table B1: Low income spending preferences",
          omit = c("country", "year"),
          covariate.labels = c("Economic orthodoxy (average, t-1)", "Economic orthodoxy (consensus, t-1)", "Economic orthodoxy (party, t-1)", "GDP (logged) (t)", "Unemployment rate (t)", "Growth (t)", "Age", "Female", "Economic orthodoxy average * consensus"),
          dep.var.labels = c("Env", "Health", "Law", "Edu", "Def", "Pension", "Unemp", "Cult", "Total"),
          add.lines=list(c("Income group", "L", "L", "L", "L", "L", "L", "L", "L", "L"),
                         c("Country-years", "60", "60", "60", "60", "60", "60","60", "60", "60", "60"),
                         c("Fixed effects", "C&Y", "C&Y", "C&Y", "C&Y", "C&Y", "C&Y", "C&Y", "C&Y", "C&Y")))

# Remove all models, only keep the dataframe
rm(list=setdiff(ls(), "df"))

# Appendix B3: spend pref mid
m1 <-   lm(spend_env ~       l_ip_per414_mean * l_ip_per414_consensus + l_ip_per414 + as.factor(country) + as.factor(year) + loggdp + unemp + growth + age + gender, weights = weight, data = subset(df, inc.group==1))
m2 <-   lm(spend_health ~    l_ip_per414_mean * l_ip_per414_consensus + l_ip_per414 + as.factor(country) + as.factor(year) + loggdp + unemp + growth + age + gender, weights = weight, data = subset(df, inc.group==1))
m3 <-   lm(spend_law ~       l_ip_per414_mean * l_ip_per414_consensus + l_ip_per414 + as.factor(country) + as.factor(year) + loggdp + unemp + growth + age + gender, weights = weight, data = subset(df, inc.group==1))
m4 <-   lm(spend_edu ~       l_ip_per414_mean * l_ip_per414_consensus + l_ip_per414 + as.factor(country) + as.factor(year) + loggdp + unemp + growth + age + gender, weights = weight, data = subset(df, inc.group==1))
m5 <-   lm(spend_def ~       l_ip_per414_mean * l_ip_per414_consensus + l_ip_per414 + as.factor(country) + as.factor(year) + loggdp + unemp + growth + age + gender, weights = weight, data = subset(df, inc.group==1))
m6 <-   lm(spend_pension ~   l_ip_per414_mean * l_ip_per414_consensus + l_ip_per414 + as.factor(country) + as.factor(year) + loggdp + unemp + growth + age + gender, weights = weight, data = subset(df, inc.group==1))
m7 <-   lm(spend_unemp ~     l_ip_per414_mean * l_ip_per414_consensus + l_ip_per414 + as.factor(country) + as.factor(year) + loggdp + unemp + growth + age + gender, weights = weight, data = subset(df, inc.group==1))
m8 <-   lm(spend_culture ~   l_ip_per414_mean * l_ip_per414_consensus + l_ip_per414 + as.factor(country) + as.factor(year) + loggdp + unemp + growth + age + gender, weights = weight, data = subset(df, inc.group==1))
m9 <-   lm(spend_total ~   l_ip_per414_mean * l_ip_per414_consensus + l_ip_per414 + as.factor(country) + as.factor(year) + loggdp + unemp + growth + age + gender, weights = weight, data = subset(df, inc.group==1))

## Robust SE: Baseline model
r1 <-  sqrt(diag(vcovHC(m1, type = "HC1")))
r2 <-  sqrt(diag(vcovHC(m2, type = "HC1")))
r3 <-  sqrt(diag(vcovHC(m3, type = "HC1")))
r4 <-  sqrt(diag(vcovHC(m4, type = "HC1")))
r5 <-  sqrt(diag(vcovHC(m5, type = "HC1")))
r6 <-  sqrt(diag(vcovHC(m6, type = "HC1")))
r7 <-  sqrt(diag(vcovHC(m7, type = "HC1")))
r8 <-  sqrt(diag(vcovHC(m8, type = "HC1")))
r9 <-  sqrt(diag(vcovHC(m9, type = "HC1")))

# Regression
stargazer(m1, m2, m3, m4, m5, m6, m7, m8, m9,
          se = list(r1, r2, r3, r4, r5, r6, r7, r8, r9),
          keep.stat=c("n", "adj.rsq"),
          type = "text", 
          #out = "../05_figures/v3/appendix.ortho.midspend.html",
          title = "Table B2: Middle income spending preferences",
          omit = c("country", "year"),
          covariate.labels = c("Economic orthodoxy (average, t-1)", "Economic orthodoxy (consensus, t-1)", "Economic orthodoxy (party, t-1)", "GDP (logged) (t)", "Unemployment rate (t)", "Growth (t)", "Age", "Female", "Economic orthodoxy average * consensus"),
          dep.var.labels = c("Env", "Health", "Law", "Edu", "Def", "Pension", "Unemp", "Cult"),
          add.lines=list(c("Income group", "M", "M", "M", "M", "M", "M", "M", "M"),
                         c("Country-years", "60", "60", "60", "60", "60", "60","60", "60", "60"),
                         c("Fixed effects", "C&Y", "C&Y", "C&Y", "C&Y", "C&Y", "C&Y", "C&Y", "C&Y", "C&Y")))

# Remove all models, only keep the dataframe
rm(list=setdiff(ls(), "df"))

#Appendix B4: spend pref high
m1 <-   lm(spend_env ~       l_ip_per414_mean * l_ip_per414_consensus + l_ip_per414 + as.factor(country) + as.factor(year) + loggdp + unemp + growth + age + gender, weights = weight, data = subset(df, inc.group==2))
m2 <-   lm(spend_health ~    l_ip_per414_mean * l_ip_per414_consensus + l_ip_per414 + as.factor(country) + as.factor(year) + loggdp + unemp + growth + age + gender, weights = weight, data = subset(df, inc.group==2))
m3 <-   lm(spend_law ~       l_ip_per414_mean * l_ip_per414_consensus + l_ip_per414 + as.factor(country) + as.factor(year) + loggdp + unemp + growth + age + gender, weights = weight, data = subset(df, inc.group==2))
m4 <-   lm(spend_edu ~       l_ip_per414_mean * l_ip_per414_consensus + l_ip_per414 + as.factor(country) + as.factor(year) + loggdp + unemp + growth + age + gender, weights = weight, data = subset(df, inc.group==2))
m5 <-   lm(spend_def ~       l_ip_per414_mean * l_ip_per414_consensus + l_ip_per414 + as.factor(country) + as.factor(year) + loggdp + unemp + growth + age + gender, weights = weight, data = subset(df, inc.group==2))
m6 <-   lm(spend_pension ~   l_ip_per414_mean * l_ip_per414_consensus + l_ip_per414 + as.factor(country) + as.factor(year) + loggdp + unemp + growth + age + gender, weights = weight, data = subset(df, inc.group==2))
m7 <-   lm(spend_unemp ~     l_ip_per414_mean * l_ip_per414_consensus + l_ip_per414 + as.factor(country) + as.factor(year) + loggdp + unemp + growth + age + gender, weights = weight, data = subset(df, inc.group==2))
m8 <-   lm(spend_culture ~   l_ip_per414_mean * l_ip_per414_consensus + l_ip_per414 + as.factor(country) + as.factor(year) + loggdp + unemp + growth + age + gender, weights = weight, data = subset(df, inc.group==2))
m9 <-   lm(spend_total ~   l_ip_per414_mean * l_ip_per414_consensus + l_ip_per414 + as.factor(country) + as.factor(year) + loggdp + unemp + growth + age + gender, weights = weight, data = subset(df, inc.group==2))

## Robust SE: Baseline model
r1 <-  sqrt(diag(vcovHC(m1, type = "HC1")))
r2 <-  sqrt(diag(vcovHC(m2, type = "HC1")))
r3 <-  sqrt(diag(vcovHC(m3, type = "HC1")))
r4 <-  sqrt(diag(vcovHC(m4, type = "HC1")))
r5 <-  sqrt(diag(vcovHC(m5, type = "HC1")))
r6 <-  sqrt(diag(vcovHC(m6, type = "HC1")))
r7 <-  sqrt(diag(vcovHC(m7, type = "HC1")))
r8 <-  sqrt(diag(vcovHC(m8, type = "HC1")))
r9 <-  sqrt(diag(vcovHC(m9, type = "HC1")))

# Regression
stargazer(m1, m2, m3, m4, m5, m6, m7, m8, m9,
          se = list(r1, r2, r3, r4, r5, r6, r7, r8, r9),
          keep.stat=c("n", "adj.rsq"),
          type = "text", 
          #out = "../05_figures/v3/appendix.ortho.highspend.html",
          title = "Table B3: High income spending preferences",
          omit = c("country", "year"),
          covariate.labels = c("Economic orthodoxy (average, t-1)", "Economic orthodoxy (consensus, t-1)", "Economic orthodoxy (party, t-1)", "GDP (logged) (t)", "Unemployment rate (t)", "Growth (t)", "Age", "Female", "Economic orthodoxy average * consensus"),
          dep.var.labels = c("Env", "Health", "Law", "Edu", "Def", "Pension", "Unemp", "Cult"),
          add.lines=list(c("Income group", "H", "H", "H", "H", "H", "H", "H", "H"),
                         c("Country-years", "60", "60", "60", "60", "60", "60","60", "60", "60"),
                         c("Fixed effects", "C&Y", "C&Y", "C&Y", "C&Y", "C&Y", "C&Y", "C&Y", "C&Y", "C&Y")))

# Remove all models, only keep the dataframe
rm(list=setdiff(ls(), "df"))

# Appendix B5: tax pref low
m1 <-   lm(tax_L ~ l_ip_per414_mean * l_ip_per414_consensus + l_ip_per414 + as.factor(country) + as.factor(year) + loggdp + unemp + growth + age + gender, weights = weight, data = subset(df, inc.group==0))
m2 <-   lm(tax_M ~ l_ip_per414_mean * l_ip_per414_consensus + l_ip_per414 + as.factor(country) + as.factor(year) + loggdp + unemp + growth + age + gender, weights = weight, data = subset(df, inc.group==0))
m3 <-   lm(tax_H ~ l_ip_per414_mean * l_ip_per414_consensus + l_ip_per414 + as.factor(country) + as.factor(year) + loggdp + unemp + growth + age + gender, weights = weight, data = subset(df, inc.group==0))
m4 <-   lm(tax_total ~ l_ip_per414_mean * l_ip_per414_consensus + l_ip_per414 + as.factor(country) + as.factor(year) + loggdp + unemp + growth + age + gender, weights = weight, data = subset(df, inc.group==0))

## Robust SE: Baseline model
r1 <-  sqrt(diag(vcovHC(m1, type = "HC1")))
r2 <-  sqrt(diag(vcovHC(m2, type = "HC1")))
r3 <-  sqrt(diag(vcovHC(m3, type = "HC1")))
r4 <-  sqrt(diag(vcovHC(m4, type = "HC1")))

# Regression
stargazer(m1, m2, m3, m4,
          se = list(r1, r2, r3, r4),
          keep.stat=c("n", "adj.rsq"),
          type = "text", 
          #out = "../05_figures/v3/appendix.ortho.lowtax.html",
          omit = c("country", "year"),
          title = "Table B4: Low income taxation preferences",
          #covariate.labels = c("Economic orthodoxy (average, t-1)", "Economic orthodoxy (consensus, t-1)", "Economic orthodoxy (party, t-1)", "GDP (logged) (t)", "Unemployment rate (t)", "Growth (t)", "Age", "Female", "Economic orthodoxy average * consensus"),
          dep.var.labels = c("Tax lower incomes", "Tax middle incomes", "Tax high incomes", "Tax total"),
          add.lines=list(c("Income group", "L", "L", "L", "L"),
                         c("Country-years", "60", "60", "60", "60"),
                         c("Fixed effects", "C&Y", "C&Y", "C&Y", "C&Y")))

# Remove all models, only keep the dataframe
rm(list=setdiff(ls(), "df"))

# Appendix B6: tax pref mid
m1 <-   lm(tax_L ~ l_ip_per414_mean * l_ip_per414_consensus + l_ip_per414 + as.factor(country) + as.factor(year) + loggdp + unemp + growth + age + gender, weights = weight, data = subset(df, inc.group==1))
m2 <-   lm(tax_M ~ l_ip_per414_mean * l_ip_per414_consensus + l_ip_per414 + as.factor(country) + as.factor(year) + loggdp + unemp + growth + age + gender, weights = weight, data = subset(df, inc.group==1))
m3 <-   lm(tax_H ~ l_ip_per414_mean * l_ip_per414_consensus + l_ip_per414 + as.factor(country) + as.factor(year) + loggdp + unemp + growth + age + gender, weights = weight, data = subset(df, inc.group==1))
m4 <-   lm(tax_total ~ l_ip_per414_mean * l_ip_per414_consensus + l_ip_per414 + as.factor(country) + as.factor(year) + loggdp + unemp + growth + age + gender, weights = weight, data = subset(df, inc.group==1))

## Robust SE: Baseline model
r1 <-  sqrt(diag(vcovHC(m1, type = "HC1")))
r2 <-  sqrt(diag(vcovHC(m2, type = "HC1")))
r3 <-  sqrt(diag(vcovHC(m3, type = "HC1")))
r4 <-  sqrt(diag(vcovHC(m4, type = "HC1")))

# Regression
stargazer(m1, m2, m3, m4,
          se = list(r1, r2, r3, r4),
          keep.stat=c("n", "adj.rsq"),
          type = "html", 
          out = "../05_figures/v3/appendix.ortho.midtax.html",
          title = "Table B5: Middle income taxation preferences",
          omit = c("country", "year"),
          covariate.labels = c("Economic orthodoxy (average, t-1)", "Economic orthodoxy (consensus, t-1)", "Economic orthodoxy (party, t-1)", "GDP (logged) (t)", "Unemployment rate (t)", "Growth (t)", "Age", "Female", "Economic orthodoxy average * consensus"),
          dep.var.labels = c("Tax lower incomes", "Tax middle incomes", "Tax high incomes", "Tax total"),
          add.lines=list(c("Income group", "M", "M", "M", "M"),
                         c("Country-years", "60", "60", "60", "60"),
                         c("Fixed effects", "C&Y", "C&Y", "C&Y", "C&Y")))

# Remove all models, only keep the dataframe
rm(list=setdiff(ls(), "df"))

# Appendix B7: tax pref high
m1 <-   lm(tax_L ~ l_ip_per414_mean * l_ip_per414_consensus + l_ip_per414 + as.factor(country) + as.factor(year) + loggdp + unemp + growth + age + gender, weights = weight, data = subset(df, inc.group==2))
m2 <-   lm(tax_M ~ l_ip_per414_mean * l_ip_per414_consensus + l_ip_per414 + as.factor(country) + as.factor(year) + loggdp + unemp + growth + age + gender, weights = weight, data = subset(df, inc.group==2))
m3 <-   lm(tax_H ~ l_ip_per414_mean * l_ip_per414_consensus + l_ip_per414 + as.factor(country) + as.factor(year) + loggdp + unemp + growth + age + gender, weights = weight, data = subset(df, inc.group==2))
m4 <-   lm(tax_total ~ l_ip_per414_mean * l_ip_per414_consensus + l_ip_per414 + as.factor(country) + as.factor(year) + loggdp + unemp + growth + age + gender, weights = weight, data = subset(df, inc.group==2))

## Robust SE: Baseline model
r1 <-  sqrt(diag(vcovHC(m1, type = "HC1")))
r2 <-  sqrt(diag(vcovHC(m2, type = "HC1")))
r3 <-  sqrt(diag(vcovHC(m3, type = "HC1")))
r4 <-  sqrt(diag(vcovHC(m4, type = "HC1")))

# Regression
stargazer(m1, m2, m3, m4,
          se = list(r1, r2, r3, r4),
          keep.stat=c("n", "adj.rsq"),
          type = "html", 
          out = "../05_figures/v3/appendix.ortho.hightax.html",
          title = "Table B6: High income taxation preferences",
          omit = c("country", "year"),
          covariate.labels = c("Economic orthodoxy (average, t-1)", "Economic orthodoxy (consensus, t-1)", "Economic orthodoxy (party, t-1)", "GDP (logged) (t)", "Unemployment rate (t)", "Growth (t)", "Age", "Female", "Economic orthodoxy average * consensus"),
          dep.var.labels = c("Tax lower incomes", "Tax middle incomes", "Tax high incomes", "Tax total"),
          add.lines=list(c("Income group", "H", "H", "H", "H"),
                         c("Country-years", "60", "60", "60", "60"),
                         c("Fixed effects", "C&Y", "C&Y", "C&Y", "C&Y")))

# Remove all models, only keep the dataframe
rm(list=setdiff(ls(), "df"))

############################################################################################################
######################## Appendix C: Results using Keynesian Demand Management #############################
############################################################################################################

# Appendix C1
m1.0 <- lm(cut_spend ~   l_ip_per409_mean * l_ip_per409_consensus + l_ip_per409 + country + year + loggdp + unemp + growth + age + gender, weights = weight, data = df)
m1 <-   lm(cut_spend ~   l_ip_per409_mean * l_ip_per409_consensus + l_ip_per409 + country + year + loggdp + unemp + growth + age + gender, weights = weight, data = subset(df, inc.group==0))
m2 <-   lm(cut_spend ~   l_ip_per409_mean * l_ip_per409_consensus + l_ip_per409 + country + year + loggdp + unemp + growth + age + gender, weights = weight, data = subset(df, inc.group==1))
m3 <-   lm(cut_spend ~   l_ip_per409_mean * l_ip_per409_consensus + l_ip_per409 + country + year + loggdp + unemp + growth + age + gender, weights = weight, data = subset(df, inc.group==2))

## Robust SE: Baseline model
r1.0 <-  sqrt(diag(vcovHC(m1.0, type = "HC1")))
r1 <-  sqrt(diag(vcovHC(m1, type = "HC1")))
r2 <-  sqrt(diag(vcovHC(m2, type = "HC1")))
r3 <-  sqrt(diag(vcovHC(m3, type = "HC1")))

# Regression
stargazer(m1.0, m1, m2, m3,
          se = list(r1.0, r1, r2, r3),
          keep.stat=c("n", "adj.rsq"),
          type = "text", 
          # out = "../05_figures/v3/appendix_c1.html",
          omit = c("country", "year"),
          covariate.labels = c("Keynesian demand management (average, t-1)", "Keynesian demand management (consensus, t-1)", "Keynesian demand management (party, t-1)", "GDP (logged) (t)", "Unemployment rate (t)", "Growth (t)", "Age", "Female", "Keynesian demand management average * consensus"),
          dep.var.labels = c("Preference for spending cuts"),
          add.lines=list(c("Income group", "All",  "L", "M", "H"),
                         c("Country-years", "60", "60", "60", "60"),
                         c("Fixed effects", "C&Y", "C&Y", "C&Y", "C&Y", "C&Y", "C&Y", "C&Y", "C&Y", "C&Y")))


# Remove all models, only keep the dataframe
rm(list=setdiff(ls(), "df"))

# Appendix C2
m1.0 <- lm(spend_polarizing ~   l_ip_per409_mean * l_ip_per409_consensus + l_ip_per409 + country + year + loggdp + unemp + growth + age + gender, weights = weight, data = df)
m1 <-   lm(spend_polarizing ~   l_ip_per409_mean * l_ip_per409_consensus + l_ip_per409 + country + year + loggdp + unemp + growth + age + gender, weights = weight, data = subset(df, inc.group==0))
m2 <-   lm(spend_polarizing ~   l_ip_per409_mean * l_ip_per409_consensus + l_ip_per409 + country + year + loggdp + unemp + growth + age + gender, weights = weight, data = subset(df, inc.group==1))
m3 <-   lm(spend_polarizing ~   l_ip_per409_mean * l_ip_per409_consensus + l_ip_per409 + country + year + loggdp + unemp + growth + age + gender, weights = weight, data = subset(df, inc.group==2))

## Robust SE: Baseline model
r1.0 <-  sqrt(diag(vcovHC(m1.0, type = "HC1")))
r1 <-  sqrt(diag(vcovHC(m1, type = "HC1")))
r2 <-  sqrt(diag(vcovHC(m2, type = "HC1")))
r3 <-  sqrt(diag(vcovHC(m3, type = "HC1")))

# Regression
stargazer(m1.0, m1, m2, m3,
          se = list(r1.0, r1, r2, r3),
          keep.stat=c("n", "adj.rsq"),
          type = "text", 
          # out = "../05_figures/v3/appendix_c2.html",
          omit = c("country", "year"),
          covariate.labels = c("Keynesian demand management (average, t-1)", "Keynesian demand management (consensus, t-1)", "Keynesian demand management (party, t-1)", "GDP (logged) (t)", "Unemployment rate (t)", "Growth (t)", "Age", "Female", "Keynesian demand management average * consensus"),
          dep.var.labels = c("Preference for spending in pensions, health and unemployment"),
          add.lines=list(c("Income group", "All",  "L", "M", "H"),
                         c("Country-years", "60", "60", "60", "60"),
                         c("Fixed effects", "C&Y", "C&Y", "C&Y", "C&Y", "C&Y", "C&Y", "C&Y", "C&Y", "C&Y")))


# Remove all models, only keep the dataframe
rm(list=setdiff(ls(), "df"))

############################################################################################################
################################## Appendix D:	Results by Left - Right ####################################
############################################################################################################

# Appendix D1: Average spending preferences across Left - Right
df %>% filter(!is.na(cut_spend) & !is.na(l_ip_per414_mean) & !is.na(l_ip_per414_consensus) & !is.na(l_ip_per414) & !is.na(country) & !is.na(year) & !is.na(unemp) & !is.na(growth) & !is.na(age) & !is.na(gender)) %>%
  group_by(lr) %>%
  summarise(mean_cut = mean(cut_spend, na.rm = T), 
            mean_spend = mean(spend_polarizing, na.rm = T)) %>%
  filter(!is.na(lr))

# Appendix D2: Linear regression models of voters’ preferences for cuts to government spending, Left - Right
m1 <-   lm(cut_spend ~   l_ip_per414_mean * l_ip_per414_consensus + l_ip_per414 + country + year + loggdp + unemp + growth + age + gender, weights = weight, data = subset(df, lr==0))
m2 <-   lm(cut_spend ~   l_ip_per414_mean * l_ip_per414_consensus + l_ip_per414 + country + year + loggdp + unemp + growth + age + gender, weights = weight, data = subset(df, lr==1))

## Robust SE: Baseline model
r1 <-  sqrt(diag(vcovHC(m1, type = "HC1")))
r2 <-  sqrt(diag(vcovHC(m2, type = "HC1")))

# Regression
stargazer(m1, m2,
          se = list(r1, r2),
          keep.stat=c("n", "adj.rsq"),
          type = "text", 
          # out = "../05_figures/v3/appendix_d1.html",
          omit = c("country", "year"),
          covariate.labels = c("Economic orthodoxy (average, t-1)", "Economic orthodoxy (consensus, t-1)", "Economic orthodoxy (party, t-1)", "GDP (logged) (t)", "Unemployment rate (t)", "Growth (t)", "Age", "Female", "Economic orthodoxy average * consensus"),
          dep.var.labels = c("Preference for spending cuts"),
          add.lines=list(c("Income group", "All",  "L", "M", "H"),
                         c("Country-years", "60", "60", "60", "60"),
                         c("Fixed effects", "C&Y", "C&Y", "C&Y", "C&Y", "C&Y", "C&Y", "C&Y", "C&Y", "C&Y")))


# Remove all models, only keep the dataframe
rm(list=setdiff(ls(), "df"))

# Appendix D3
m1 <-   lm(spend_polarizing ~   l_ip_per414_mean * l_ip_per414_consensus + l_ip_per414 + country + year + loggdp + unemp + growth + age + gender, weights = weight, data = subset(df, lr==0))
m2 <-   lm(spend_polarizing ~   l_ip_per414_mean * l_ip_per414_consensus + l_ip_per414 + country + year + loggdp + unemp + growth + age + gender, weights = weight, data = subset(df, lr==1))

## Robust SE: Baseline model
r1 <-  sqrt(diag(vcovHC(m1, type = "HC1")))
r2 <-  sqrt(diag(vcovHC(m2, type = "HC1")))

# Regression
stargazer(m1, m2,
          se = list(r1, r2),
          keep.stat=c("n", "adj.rsq"),
          type = "text", 
          # out = "../05_figures/v3/appendix_d2.html",
          omit = c("country", "year"),
          covariate.labels = c("Economic orthodoxy (average, t-1)", "Economic orthodoxy (consensus, t-1)", "Economic orthodoxy (party, t-1)", "GDP (logged) (t)", "Unemployment rate (t)", "Growth (t)", "Age", "Female", "Economic orthodoxy average * consensus"),
          dep.var.labels = c("Preference for spending in pensions, health and unemployment"),
          add.lines=list(c("Income group", "All",  "L", "M", "H"),
                         c("Country-years", "60", "60", "60", "60"),
                         c("Fixed effects", "C&Y", "C&Y", "C&Y", "C&Y", "C&Y", "C&Y", "C&Y", "C&Y", "C&Y")))


# Remove all models, only keep the dataframe
rm(list=setdiff(ls(), "df"))

############################################################################################################
####################### Appendix E: Results by education and information groups ############################
############################################################################################################

# Appendix E1: Average spending preferences across Left - Right
df %>% filter(!is.na(cut_spend) & !is.na(l_ip_per414_mean) & !is.na(l_ip_per414_consensus) & !is.na(l_ip_per414) & !is.na(country) & !is.na(year) & !is.na(unemp) & !is.na(growth) & !is.na(age) & !is.na(gender)) %>%
  mutate(educat = case_when(educ_level %in% c(0,1,2) ~ 1,
                            educ_level %in% c(3) ~ 2,
                            educ_level %in% c(4,5,6) ~ 3)) %>%
  group_by(educat) %>%
  summarise(mean_cut = mean(cut_spend, na.rm = T))

df %>% filter(!is.na(cut_spend) & !is.na(l_ip_per414_mean) & !is.na(l_ip_per414_consensus) & !is.na(l_ip_per414) & !is.na(country) & !is.na(year) & !is.na(unemp) & !is.na(growth) & !is.na(age) & !is.na(gender)) %>%
  mutate(informedcat = case_when(informed %in% c(1,2) ~ 1,
                                 informed %in% c(3) ~ 2,
                                 informed %in% c(4,5) ~ 3)) %>%
  group_by(informedcat) %>%
  summarise(mean_cut = mean(cut_spend, na.rm = T))

# Appendix E2: Economic orthodoxy Other categories
# Education
m1 <-   lm(cut_spend ~   l_ip_per414_mean * l_ip_per414_consensus + l_ip_per414 + country + year + loggdp + unemp + growth + age + gender, weights = weight, data = subset(df, educ_level %in% c(0,1,2)))
m2 <-   lm(cut_spend ~   l_ip_per414_mean * l_ip_per414_consensus + l_ip_per414 + country + year + loggdp + unemp + growth + age + gender, weights = weight, data = subset(df, educ_level %in% c(3)))
m3 <-   lm(cut_spend ~   l_ip_per414_mean * l_ip_per414_consensus + l_ip_per414 + country + year + loggdp + unemp + growth + age + gender, weights = weight, data = subset(df, educ_level %in% c(4,5,6)))

# Informed
m4 <-   lm(cut_spend ~   l_ip_per414_mean * l_ip_per414_consensus + l_ip_per414 + country + year + loggdp + unemp + growth + age + gender, weights = weight, data = subset(df, informed %in% c(1,2)))
m5 <-   lm(cut_spend ~   l_ip_per414_mean * l_ip_per414_consensus + l_ip_per414 + country + year + loggdp + unemp + growth + age + gender, weights = weight, data = subset(df, informed %in% c(3)))
m6 <-   lm(cut_spend ~   l_ip_per414_mean * l_ip_per414_consensus + l_ip_per414 + country + year + loggdp + unemp + growth + age + gender, weights = weight, data = subset(df, informed %in% c(4,5)))

## Robust SE: Baseline model
r1 <-  sqrt(diag(vcovHC(m1, type = "HC1")))
r2 <-  sqrt(diag(vcovHC(m2, type = "HC1")))
r3 <-  sqrt(diag(vcovHC(m3, type = "HC1")))
r4 <-  sqrt(diag(vcovHC(m4, type = "HC1")))
r5 <-  sqrt(diag(vcovHC(m5, type = "HC1")))
r6 <-  sqrt(diag(vcovHC(m6, type = "HC1")))

# Regression
stargazer(m1, m2, m3, m4, m5, m6,
          se = list(r1, r2, r3, r4, r5, r6),
          keep.stat=c("n", "adj.rsq"),
          type = "text", 
          #out = "../05_figures/v3/ortho_other_categories.html",
          omit = c("country", "year"),
          covariate.labels = c("Economic orthodoxy (average, t-1)", "Economic orthodoxy (consensus, t-1)", "Economic orthodoxy (party, t-1)", "GDP (logged) (t)", "Unemployment rate (t)", "Growth (t)", "Age", "Female", "Economic orthodoxy average * consensus"),
          dep.var.labels = c("Preference for spending cuts"),
          add.lines=list(c("Group", "Low Edu", "Mid Edu",  "High Edu", "Low Informed", "Mid Informed", "High Informed"),
                         c("Country-years", "60", "60", "60", "60"),
                         c("Fixed effects", "C&Y", "C&Y", "C&Y", "C&Y", "C&Y", "C&Y", "C&Y", "C&Y", "C&Y")))

# Remove all models, only keep the dataframe
rm(list=setdiff(ls(), "df"))

# Appendix E3: Other categories
# Education
m1 <-   lm(spend_polarizing ~   l_ip_per414_mean * l_ip_per414_consensus + l_ip_per414 + country + year + loggdp + unemp + growth + age + gender, weights = weight, data = subset(df, educ_level %in% c(0,1,2)))
m2 <-   lm(spend_polarizing ~   l_ip_per414_mean * l_ip_per414_consensus + l_ip_per414 + country + year + loggdp + unemp + growth + age + gender, weights = weight, data = subset(df, educ_level %in% c(3)))
m3 <-   lm(spend_polarizing ~   l_ip_per414_mean * l_ip_per414_consensus + l_ip_per414 + country + year + loggdp + unemp + growth + age + gender, weights = weight, data = subset(df, educ_level %in% c(4,5,6)))

# Informed
m4 <-   lm(spend_polarizing ~   l_ip_per414_mean * l_ip_per414_consensus + l_ip_per414 + country + year + loggdp + unemp + growth + age + gender, weights = weight, data = subset(df, informed %in% c(1,2)))
m5 <-   lm(spend_polarizing ~   l_ip_per414_mean * l_ip_per414_consensus + l_ip_per414 + country + year + loggdp + unemp + growth + age + gender, weights = weight, data = subset(df, informed %in% c(3)))
m6 <-   lm(spend_polarizing ~   l_ip_per414_mean * l_ip_per414_consensus + l_ip_per414 + country + year + loggdp + unemp + growth + age + gender, weights = weight, data = subset(df, informed %in% c(4,5)))

## Robust SE: Baseline model
r1 <-  sqrt(diag(vcovHC(m1, type = "HC1")))
r2 <-  sqrt(diag(vcovHC(m2, type = "HC1")))
r3 <-  sqrt(diag(vcovHC(m3, type = "HC1")))
r4 <-  sqrt(diag(vcovHC(m4, type = "HC1")))
r5 <-  sqrt(diag(vcovHC(m5, type = "HC1")))
r6 <-  sqrt(diag(vcovHC(m6, type = "HC1")))

# Regression
stargazer(m1, m2, m3, m4, m5, m6,
          se = list(r1, r2, r3, r4, r5, r6),
          keep.stat=c("n", "adj.rsq"),
          type = "text", 
          #out = "../05_figures/v3/keynes_other_categories.html",
          omit = c("country", "year"),
          covariate.labels = c("Keynesian demand management (average, t-1)", "Keynesian demand management (consensus, t-1)", "Keynesian demand management (party, t-1)", "GDP (logged) (t)", "Unemployment rate (t)", "Growth (t)", "Age", "Female", "Keynesian demand management average * consensus"),
          dep.var.labels = c("Preference for spending cuts"),
          add.lines=list(c("Group", "Low Edu", "Mid Edu",  "High Edu", "Low Informed", "Mid Informed", "High Informed"),
                         c("Country-years", "60", "60", "60", "60"),
                         c("Fixed effects", "C&Y", "C&Y", "C&Y", "C&Y", "C&Y", "C&Y", "C&Y", "C&Y", "C&Y")))


# Remove all models, only keep the dataframe
rm(list=setdiff(ls(), "df"))

############################################################################################################
####################################### Appendix F: Robustness tests #######################################
############################################################################################################

#F1 Linear regression models of voters’ preferences for cuts to government spending, subset survey within two years of election
m1.0 <- lm(cut_spend ~   l_ip_per414_mean * l_ip_per414_consensus + l_ip_per414 + country + year + loggdp + unemp + growth + age + gender, weights = weight, data = subset(df, yeardistance<=2))
m1 <-   lm(cut_spend ~   l_ip_per414_mean * l_ip_per414_consensus + l_ip_per414 + country + year + loggdp + unemp + growth + age + gender, weights = weight, data = subset(df, inc.group==0 & yeardistance<=2))
m2 <-   lm(cut_spend ~   l_ip_per414_mean * l_ip_per414_consensus + l_ip_per414 + country + year + loggdp + unemp + growth + age + gender, weights = weight, data = subset(df, inc.group==1 & yeardistance<=2))
m3 <-   lm(cut_spend ~   l_ip_per414_mean * l_ip_per414_consensus + l_ip_per414 + country + year + loggdp + unemp + growth + age + gender, weights = weight, data = subset(df, inc.group==2 & yeardistance<=2))

## Robust SE: Baseline model
r1.0 <-  sqrt(diag(vcovHC(m1.0, type = "HC1")))
r1 <-  sqrt(diag(vcovHC(m1, type = "HC1")))
r2 <-  sqrt(diag(vcovHC(m2, type = "HC1")))
r3 <-  sqrt(diag(vcovHC(m3, type = "HC1")))

# Regression
stargazer(m1.0, m1, m2, m3,
          se = list(r1.0, r1, r2, r3),
          keep.stat=c("n", "adj.rsq"),
          type = "text", 
          # out = "../05_figures/v3/main1.html",
          omit = c("country", "year"),
          covariate.labels = c("Economic orthodoxy (average, t-1)", "Economic orthodoxy (consensus, t-1)", "Economic orthodoxy (party, t-1)", "GDP (logged) (t)", "Unemployment rate (t)", "Growth (t)", "Age", "Female", "Economic orthodoxy average * consensus"),
          dep.var.labels = c("Preference for spending cuts"),
          add.lines=list(c("Income group", "All",  "L", "M", "H"),
                         c("Country-years", "60", "60", "60", "60"),
                         c("Fixed effects", "C&Y", "C&Y", "C&Y", "C&Y", "C&Y", "C&Y", "C&Y", "C&Y", "C&Y")))

# Remove all models, only keep the dataframe
rm(list=setdiff(ls(), "df"))

#F2: Linear regression models of voters’ preferences for spending in pensions, health and unemployment, subset survey within two years of election
m1.0 <- lm(spend_polarizing ~   l_ip_per414_mean * l_ip_per414_consensus + l_ip_per414 + country + year + loggdp + unemp + growth + age + gender, weights = weight, data = subset(df, yeardistance<=2))
m1 <-   lm(spend_polarizing ~   l_ip_per414_mean * l_ip_per414_consensus + l_ip_per414 + country + year + loggdp + unemp + growth + age + gender, weights = weight, data = subset(df, inc.group==0 & yeardistance<=2))
m2 <-   lm(spend_polarizing ~   l_ip_per414_mean * l_ip_per414_consensus + l_ip_per414 + country + year + loggdp + unemp + growth + age + gender, weights = weight, data = subset(df, inc.group==1 & yeardistance<=2))
m3 <-   lm(spend_polarizing ~   l_ip_per414_mean * l_ip_per414_consensus + l_ip_per414 + country + year + loggdp + unemp + growth + age + gender, weights = weight, data = subset(df, inc.group==2 & yeardistance<=2))

## Robust SE: Baseline model
r1.0 <-  sqrt(diag(vcovHC(m1.0, type = "HC1")))
r1 <-  sqrt(diag(vcovHC(m1, type = "HC1")))
r2 <-  sqrt(diag(vcovHC(m2, type = "HC1")))
r3 <-  sqrt(diag(vcovHC(m3, type = "HC1")))

# Regression
stargazer(m1.0, m1, m2, m3,
          se = list(r1.0, r1, r2, r3),
          keep.stat=c("n", "adj.rsq"),
          type = "html", 
          #out = "f2.html",
          omit = c("country", "year"),
          covariate.labels = c("Economic orthodoxy (average, t-1)", "Economic orthodoxy (consensus, t-1)", "Economic orthodoxy (party, t-1)", "GDP (logged) (t)", "Unemployment rate (t)", "Growth (t)", "Age", "Female", "Economic orthodoxy average * consensus"),
          dep.var.labels = c("Preference for spending in pensions, health and unemployment"),
          add.lines=list(c("Income group", "All",  "L", "M", "H"),
                         c("Country-years", "60", "60", "60", "60"),
                         c("Fixed effects", "C&Y", "C&Y", "C&Y", "C&Y", "C&Y", "C&Y", "C&Y", "C&Y", "C&Y")))


# F4: Linear regression models of voters’ preferences for cuts to government spending, austerity composite
m1.0 <- lm(cut_spend ~   l_ip_perall_mean * l_ip_perall_consensus + l_ip_perall + country + year + loggdp + unemp + growth + age + gender, weights = weight, data = df)
m1 <-   lm(cut_spend ~   l_ip_perall_mean * l_ip_perall_consensus + l_ip_perall + country + year + loggdp + unemp + growth + age + gender, weights = weight, data = subset(df, inc.group==0))
m2 <-   lm(cut_spend ~   l_ip_perall_mean * l_ip_perall_consensus + l_ip_perall + country + year + loggdp + unemp + growth + age + gender, weights = weight, data = subset(df, inc.group==1))
m3 <-   lm(cut_spend ~   l_ip_perall_mean * l_ip_perall_consensus + l_ip_perall + country + year + loggdp + unemp + growth + age + gender, weights = weight, data = subset(df, inc.group==2))

## Robust SE: Baseline model
r1.0 <-  sqrt(diag(vcovHC(m1.0, type = "HC1")))
r1 <-  sqrt(diag(vcovHC(m1, type = "HC1")))
r2 <-  sqrt(diag(vcovHC(m2, type = "HC1")))
r3 <-  sqrt(diag(vcovHC(m3, type = "HC1")))

# Regression
stargazer(m1.0, m1, m2, m3,
          se = list(r1.0, r1, r2, r3),
          keep.stat=c("n", "adj.rsq"),
          type = "text", 
          # out = "../05_figures/v3/main1.html",
          omit = c("country", "year"),
          covariate.labels = c("Austerity composite (average, t-1)", "Austerity composite (consensus, t-1)", "Austerity composite (party, t-1)", "GDP (logged) (t)", "Unemployment rate (t)", "Growth (t)", "Age", "Female", "Economic orthodoxy average * consensus"),
          dep.var.labels = c("Preference for spending cuts"),
          add.lines=list(c("Income group", "All",  "L", "M", "H"),
                         c("Country-years", "60", "60", "60", "60"),
                         c("Fixed effects", "C&Y", "C&Y", "C&Y", "C&Y", "C&Y", "C&Y", "C&Y", "C&Y", "C&Y")))


# Remove all models, only keep the dataframe
rm(list=setdiff(ls(), "df"))

# F5: Linear regression models of voters’ preferences for spending in pensions, health and unemployment, austerity composite
m1.0 <- lm(spend_polarizing ~   l_ip_perall_mean * l_ip_perall_consensus + l_ip_perall + country + year + loggdp + unemp + growth + age + gender, weights = weight, data = df)
m1 <-   lm(spend_polarizing ~   l_ip_perall_mean * l_ip_perall_consensus + l_ip_perall + country + year + loggdp + unemp + growth + age + gender, weights = weight, data = subset(df, inc.group==0))
m2 <-   lm(spend_polarizing ~   l_ip_perall_mean * l_ip_perall_consensus + l_ip_perall + country + year + loggdp + unemp + growth + age + gender, weights = weight, data = subset(df, inc.group==1))
m3 <-   lm(spend_polarizing ~   l_ip_perall_mean * l_ip_perall_consensus + l_ip_perall + country + year + loggdp + unemp + growth + age + gender, weights = weight, data = subset(df, inc.group==2))

## Robust SE: Baseline model
r1.0 <-  sqrt(diag(vcovHC(m1.0, type = "HC1")))
r1 <-  sqrt(diag(vcovHC(m1, type = "HC1")))
r2 <-  sqrt(diag(vcovHC(m2, type = "HC1")))
r3 <-  sqrt(diag(vcovHC(m3, type = "HC1")))

# Regression
stargazer(m1.0, m1, m2, m3,
          se = list(r1.0, r1, r2, r3),
          keep.stat=c("n", "adj.rsq"),
          type = "text", 
          # out = "../05_figures/v3/main1.html",
          omit = c("country", "year"),
          covariate.labels = c("Austerity composite (average, t-1)", "Austerity composite (consensus, t-1)", "Austerity composite (party, t-1)", "GDP (logged) (t)", "Unemployment rate (t)", "Growth (t)", "Age", "Female", "Economic orthodoxy average * consensus"),
          dep.var.labels = c("Preference for spending cuts"),
          add.lines=list(c("Income group", "All",  "L", "M", "H"),
                         c("Country-years", "60", "60", "60", "60"),
                         c("Fixed effects", "C&Y", "C&Y", "C&Y", "C&Y", "C&Y", "C&Y", "C&Y", "C&Y", "C&Y")))


# Remove all models, only keep the dataframe
rm(list=setdiff(ls(), "df"))

# F6: Linear regression models of voters’ preferences for cuts to government spending, change in party position
m1.0 <- lm(cut_spend ~   per414_mean_change * per414_consensus_change + ip_per414_change + country + year + loggdp + unemp + growth + age + gender, weights = weight, data = df)
m1 <-   lm(cut_spend ~   per414_mean_change * per414_consensus_change + ip_per414_change + country + year + loggdp + unemp + growth + age + gender, weights = weight, data = subset(df, inc.group==0))
m2 <-   lm(cut_spend ~   per414_mean_change * per414_consensus_change + ip_per414_change + country + year + loggdp + unemp + growth + age + gender, weights = weight, data = subset(df, inc.group==1))
m3 <-   lm(cut_spend ~   per414_mean_change * per414_consensus_change + ip_per414_change + country + year + loggdp + unemp + growth + age + gender, weights = weight, data = subset(df, inc.group==2))

## Robust SE: Baseline model
r1.0 <-  sqrt(diag(vcovHC(m1.0, type = "HC1")))
r1 <-  sqrt(diag(vcovHC(m1, type = "HC1")))
r2 <-  sqrt(diag(vcovHC(m2, type = "HC1")))
r3 <-  sqrt(diag(vcovHC(m3, type = "HC1")))

# Regression
stargazer(m1.0, m1, m2, m3,
          se = list(r1.0, r1, r2, r3),
          keep.stat=c("n", "adj.rsq"),
          type = "text", 
          # out = "../05_figures/v3/main1.html",
          omit = c("country", "year"),
          covariate.labels = c("Economic orthodoxy (average, change)", "Economic orthodoxy (consensus, change)", "Economic orthodoxy (party, change)", "GDP (logged) (t)", "Unemployment rate (t)", "Growth (t)", "Age", "Female", "Economic orthodoxy average * consensus"),
          dep.var.labels = c("Preference for spending in pensions, health and unemployment"),
          add.lines=list(c("Income group", "All",  "L", "M", "H"),
                         c("Country-years", "60", "60", "60", "60"),
                         c("Fixed effects", "C&Y", "C&Y", "C&Y", "C&Y", "C&Y", "C&Y", "C&Y", "C&Y", "C&Y")))


# Remove all models, only keep the dataframe
rm(list=setdiff(ls(), "df"))

# F7: Linear regression models of voters’ preferences for spending in pensions, health and unemployment, change in party position
m1.0 <- lm(spend_polarizing ~   per414_mean_change * per414_consensus_change + ip_per414_change + country + year + loggdp + unemp + growth + age + gender, weights = weight, data = df)
m1 <-   lm(spend_polarizing ~   per414_mean_change * per414_consensus_change + ip_per414_change + country + year + loggdp + unemp + growth + age + gender, weights = weight, data = subset(df, inc.group==0))
m2 <-   lm(spend_polarizing ~   per414_mean_change * per414_consensus_change + ip_per414_change + country + year + loggdp + unemp + growth + age + gender, weights = weight, data = subset(df, inc.group==1))
m3 <-   lm(spend_polarizing ~   per414_mean_change * per414_consensus_change + ip_per414_change + country + year + loggdp + unemp + growth + age + gender, weights = weight, data = subset(df, inc.group==2))

## Robust SE: Baseline model
r1.0 <-  sqrt(diag(vcovHC(m1.0, type = "HC1")))
r1 <-  sqrt(diag(vcovHC(m1, type = "HC1")))
r2 <-  sqrt(diag(vcovHC(m2, type = "HC1")))
r3 <-  sqrt(diag(vcovHC(m3, type = "HC1")))

# Regression
stargazer(m1.0, m1, m2, m3,
          se = list(r1.0, r1, r2, r3),
          keep.stat=c("n", "adj.rsq"),
          type = "text", 
          # out = "../05_figures/v3/main1.html",
          omit = c("country", "year"),
          covariate.labels = c("Economic orthodoxy (average, change)", "Economic orthodoxy (consensus, change)", "Economic orthodoxy (party, change)", "GDP (logged) (t)", "Unemployment rate (t)", "Growth (t)", "Age", "Female", "Economic orthodoxy average * consensus"),
          dep.var.labels = c("Preference for spending in pensions, health and unemployment"),
          add.lines=list(c("Income group", "All",  "L", "M", "H"),
                         c("Country-years", "60", "60", "60", "60"),
                         c("Fixed effects", "C&Y", "C&Y", "C&Y", "C&Y", "C&Y", "C&Y", "C&Y", "C&Y", "C&Y")))


# Remove all models, only keep the dataframe
rm(list=setdiff(ls(), "df"))

# F8: Linear regression models of voters’ preferences for cuts to government spending, different lags and leads
m1 <-   lm(cut_spend ~   ip_per414_mean * ip_per414_consensus + ip_per414 + country + year + loggdp + unemp + growth + age + gender, weights = weight, data = df)
m2 <-   lm(cut_spend ~   ip_per414_mean * ip_per414_consensus + ip_per414 + country + year + loggdp + unemp + growth + age + gender, weights = weight, data = subset(df, inc.group==0))
m3 <-   lm(cut_spend ~   ip_per414_mean * ip_per414_consensus + ip_per414 + country + year + loggdp + unemp + growth + age + gender, weights = weight, data = subset(df, inc.group==1))
m4 <-   lm(cut_spend ~   ip_per414_mean * ip_per414_consensus + ip_per414 + country + year + loggdp + unemp + growth + age + gender, weights = weight, data = subset(df, inc.group==2))
m5 <-   lm(cut_spend ~   l2_ip_per414_mean * l2_ip_per414_consensus + l2_ip_per414 + country + year + loggdp + unemp + growth + age + gender, weights = weight, data = df)
m6 <-   lm(cut_spend ~   l2_ip_per414_mean * l2_ip_per414_consensus + l2_ip_per414 + country + year + loggdp + unemp + growth + age + gender, weights = weight, data = subset(df, inc.group==0))
m7 <-   lm(cut_spend ~   l2_ip_per414_mean * l2_ip_per414_consensus + l2_ip_per414 + country + year + loggdp + unemp + growth + age + gender, weights = weight, data = subset(df, inc.group==1))
m8 <-   lm(cut_spend ~   l2_ip_per414_mean * l2_ip_per414_consensus + l2_ip_per414 + country + year + loggdp + unemp + growth + age + gender, weights = weight, data = subset(df, inc.group==2))

## Robust SE: Baseline model
r1 <-  sqrt(diag(vcovHC(m1, type = "HC1")))
r2 <-  sqrt(diag(vcovHC(m2, type = "HC1")))
r3 <-  sqrt(diag(vcovHC(m3, type = "HC1")))
r4 <-  sqrt(diag(vcovHC(m4, type = "HC1")))
r5 <-  sqrt(diag(vcovHC(m5, type = "HC1")))
r6 <-  sqrt(diag(vcovHC(m6, type = "HC1")))
r7 <-  sqrt(diag(vcovHC(m7, type = "HC1")))
r8 <-  sqrt(diag(vcovHC(m8, type = "HC1")))

# Regression
stargazer(m1, m2, m3, m4, m5, m6, m7, m8,
          se = list(r1, r2, r3, r4, r5, r6, r7, r8),
          keep.stat=c("n", "adj.rsq"),
          type = "text", 
          # out = "../05_figures/v3/f9.html",
          omit = c("country", "year"),
          covariate.labels = c("Economic orthodoxy (average, t)", "Economic orthodoxy (consensus, t)", "Economic orthodoxy (party, t)", "Economic orthodoxy (average, t-2)", "Economic orthodoxy (consensus, t-2)", "Economic orthodoxy (party, t-2)", "GDP (logged) (t)", "Unemployment rate (t)", "Growth (t)", "Age", "Female", "Economic orthodoxy average (t) * consensus (t)", "Economic orthodoxy average (t-2) * consensus (t-2)"),
          dep.var.labels = c("Preference for spending cuts"),
          add.lines=list(c("Income group", "All",  "L", "M", "H"),
                         c("Country-years", "60", "60", "60", "60"),
                         c("Fixed effects", "C&Y", "C&Y", "C&Y", "C&Y", "C&Y", "C&Y", "C&Y", "C&Y", "C&Y")))

# F9: Linear regression models of voters’ preferences for spending in pensions, health and unemployment, different lags and leads
m1 <-   lm(spend_polarizing ~   ip_per414_mean * ip_per414_consensus + ip_per414 + country + year + loggdp + unemp + growth + age + gender, weights = weight, data = df)
m2 <-   lm(spend_polarizing ~   ip_per414_mean * ip_per414_consensus + ip_per414 + country + year + loggdp + unemp + growth + age + gender, weights = weight, data = subset(df, inc.group==0))
m3 <-   lm(spend_polarizing ~   ip_per414_mean * ip_per414_consensus + ip_per414 + country + year + loggdp + unemp + growth + age + gender, weights = weight, data = subset(df, inc.group==1))
m4 <-   lm(spend_polarizing ~   ip_per414_mean * ip_per414_consensus + ip_per414 + country + year + loggdp + unemp + growth + age + gender, weights = weight, data = subset(df, inc.group==2))
m5 <-   lm(spend_polarizing ~   l2_ip_per414_mean * l2_ip_per414_consensus + l2_ip_per414 + country + year + loggdp + unemp + growth + age + gender, weights = weight, data = df)
m6 <-   lm(spend_polarizing ~   l2_ip_per414_mean * l2_ip_per414_consensus + l2_ip_per414 + country + year + loggdp + unemp + growth + age + gender, weights = weight, data = subset(df, inc.group==0))
m7 <-   lm(spend_polarizing ~   l2_ip_per414_mean * l2_ip_per414_consensus + l2_ip_per414 + country + year + loggdp + unemp + growth + age + gender, weights = weight, data = subset(df, inc.group==1))
m8 <-   lm(spend_polarizing ~   l2_ip_per414_mean * l2_ip_per414_consensus + l2_ip_per414 + country + year + loggdp + unemp + growth + age + gender, weights = weight, data = subset(df, inc.group==2))

## Robust SE: Baseline model
r1 <-  sqrt(diag(vcovHC(m1, type = "HC1")))
r2 <-  sqrt(diag(vcovHC(m2, type = "HC1")))
r3 <-  sqrt(diag(vcovHC(m3, type = "HC1")))
r4 <-  sqrt(diag(vcovHC(m4, type = "HC1")))
r5 <-  sqrt(diag(vcovHC(m5, type = "HC1")))
r6 <-  sqrt(diag(vcovHC(m6, type = "HC1")))
r7 <-  sqrt(diag(vcovHC(m7, type = "HC1")))
r8 <-  sqrt(diag(vcovHC(m8, type = "HC1")))

# Regression
stargazer(m1, m2, m3, m4, m5, m6, m7, m8,
          se = list(r1, r2, r3, r4, r5, r6, r7, r8),
          keep.stat=c("n", "adj.rsq"),
          type = "text", 
          # out = "../05_figures/v3/f9.html",
          omit = c("country", "year"),
          covariate.labels = c("Economic orthodoxy (average, t)", "Economic orthodoxy (consensus, t)", "Economic orthodoxy (party, t)", "Economic orthodoxy (average, t-2)", "Economic orthodoxy (consensus, t-2)", "Economic orthodoxy (party, t-2)", "GDP (logged) (t)", "Unemployment rate (t)", "Growth (t)", "Age", "Female", "Economic orthodoxy average (t) * consensus (t)", "Economic orthodoxy average (t-2) * consensus (t-2)"),
          dep.var.labels = c("Preference for spending in pensions, health and unemployment"),
          add.lines=list(c("Income group", "All",  "L", "M", "H"),
                         c("Country-years", "60", "60", "60", "60"),
                         c("Fixed effects", "C&Y", "C&Y", "C&Y", "C&Y", "C&Y", "C&Y", "C&Y", "C&Y", "C&Y")))


# F10: Linear regression models of voters’ preferences for cuts to government spending, excluding country and year dummies
m1 <-   lm(cut_spend ~   l_ip_per414_mean * l_ip_per414_consensus + l_ip_per414 + country + loggdp + unemp + growth + age + gender, weights = weight, data = df)
m2 <-   lm(cut_spend ~   l_ip_per414_mean * l_ip_per414_consensus + l_ip_per414 + country + loggdp + unemp + growth + age + gender, weights = weight, data = subset(df, inc.group==0))
m3 <-   lm(cut_spend ~   l_ip_per414_mean * l_ip_per414_consensus + l_ip_per414 + country + loggdp + unemp + growth + age + gender, weights = weight, data = subset(df, inc.group==1))
m4 <-   lm(cut_spend ~   l_ip_per414_mean * l_ip_per414_consensus + l_ip_per414 + country + loggdp + unemp + growth + age + gender, weights = weight, data = subset(df, inc.group==2))
m5 <-   lm(cut_spend ~   l_ip_per414_mean * l_ip_per414_consensus + l_ip_per414 + year + loggdp + unemp + growth + age + gender, weights = weight, data = df)
m6 <-   lm(cut_spend ~   l_ip_per414_mean * l_ip_per414_consensus + l_ip_per414 + year + loggdp + unemp + growth + age + gender, weights = weight, data = subset(df, inc.group==0))
m7 <-   lm(cut_spend ~   l_ip_per414_mean * l_ip_per414_consensus + l_ip_per414 + year + loggdp + unemp + growth + age + gender, weights = weight, data = subset(df, inc.group==1))
m8 <-   lm(cut_spend ~   l_ip_per414_mean * l_ip_per414_consensus + l_ip_per414 + year + loggdp + unemp + growth + age + gender, weights = weight, data = subset(df, inc.group==2))


## Robust SE: Baseline model
r1 <-  sqrt(diag(vcovHC(m1, type = "HC1")))
r2 <-  sqrt(diag(vcovHC(m2, type = "HC1")))
r3 <-  sqrt(diag(vcovHC(m3, type = "HC1")))
r4 <-  sqrt(diag(vcovHC(m4, type = "HC1")))
r5 <-  sqrt(diag(vcovHC(m5, type = "HC1")))
r6 <-  sqrt(diag(vcovHC(m6, type = "HC1")))
r7 <-  sqrt(diag(vcovHC(m7, type = "HC1")))
r8 <-  sqrt(diag(vcovHC(m8, type = "HC1")))


# Regression
stargazer(m1, m2, m3, m4, m5, m6, m7, m8,
          se = list(r1, r2, r3, r4, r5, r6, r7, r8),
          keep.stat=c("n", "adj.rsq"),
          type = "text", 
          # out = "../05_figures/v3/appendix_f10.html",
          omit = c("country", "year"),
          covariate.labels = c("Economic orthodoxy (average, t-1)", "Economic orthodoxy (consensus, t-1)", "GDP (logged) (t)", "Unemployment rate (t)", "Growth (t)", "Age", "Female", "Economic orthodoxy average * consensus"),
          dep.var.labels = c("Preference for spending cuts"),
          add.lines=list(c("Income group", "All",  "L", "M", "H", "All",  "L", "M", "H"),
                         c("Country-years", "60", "60", "60", "60"),
                         c("Fixed effects", "C", "C", "C", "C", "Y", "Y", "Y", "Y", "Y")))

# F11: Linear regression models of voters’ preferences for spending in pensions, health and unemployment, excluding country and year dummies
m1 <-   lm(spend_polarizing ~   l_ip_per414_mean * l_ip_per414_consensus + l_ip_per414 + country + loggdp + unemp + growth + age + gender, weights = weight, data = df)
m2 <-   lm(spend_polarizing ~   l_ip_per414_mean * l_ip_per414_consensus + l_ip_per414 + country + loggdp + unemp + growth + age + gender, weights = weight, data = subset(df, inc.group==0))
m3 <-   lm(spend_polarizing ~   l_ip_per414_mean * l_ip_per414_consensus + l_ip_per414 + country + loggdp + unemp + growth + age + gender, weights = weight, data = subset(df, inc.group==1))
m4 <-   lm(spend_polarizing ~   l_ip_per414_mean * l_ip_per414_consensus + l_ip_per414 + country + loggdp + unemp + growth + age + gender, weights = weight, data = subset(df, inc.group==2))
m5 <-   lm(spend_polarizing ~   l_ip_per414_mean * l_ip_per414_consensus + l_ip_per414 + year + loggdp + unemp + growth + age + gender, weights = weight, data = df)
m6 <-   lm(spend_polarizing ~   l_ip_per414_mean * l_ip_per414_consensus + l_ip_per414 + year + loggdp + unemp + growth + age + gender, weights = weight, data = subset(df, inc.group==0))
m7 <-   lm(spend_polarizing ~   l_ip_per414_mean * l_ip_per414_consensus + l_ip_per414 + year + loggdp + unemp + growth + age + gender, weights = weight, data = subset(df, inc.group==1))
m8 <-   lm(spend_polarizing ~   l_ip_per414_mean * l_ip_per414_consensus + l_ip_per414 + year + loggdp + unemp + growth + age + gender, weights = weight, data = subset(df, inc.group==2))

## Robust SE: Baseline model
r1 <-  sqrt(diag(vcovHC(m1, type = "HC1")))
r2 <-  sqrt(diag(vcovHC(m2, type = "HC1")))
r3 <-  sqrt(diag(vcovHC(m3, type = "HC1")))
r4 <-  sqrt(diag(vcovHC(m4, type = "HC1")))
r5 <-  sqrt(diag(vcovHC(m5, type = "HC1")))
r6 <-  sqrt(diag(vcovHC(m6, type = "HC1")))
r7 <-  sqrt(diag(vcovHC(m7, type = "HC1")))
r8 <-  sqrt(diag(vcovHC(m8, type = "HC1")))


# Regression
stargazer(m1, m2, m3, m4, m5, m6, m7, m8,
          se = list(r1, r2, r3, r4, r5, r6, r7, r8),
          keep.stat=c("n", "adj.rsq"),
          type = "text", 
          # out = "../05_figures/v3/appendix_f11.html",
          omit = c("country", "year"),
          covariate.labels = c("Economic orthodoxy (average, t-1)", "Economic orthodoxy (consensus, t-1)", "GDP (logged) (t)", "Unemployment rate (t)", "Growth (t)", "Age", "Female", "Economic orthodoxy average * consensus"),
          dep.var.labels = c("Preference for spending in pensions, health and unemployment"),
          add.lines=list(c("Income group", "All",  "L", "M", "H", "All",  "L", "M", "H"),
                         c("Country-years", "60", "60", "60", "60"),
                         c("Fixed effects", "C", "C", "C", "C", "Y", "Y", "Y", "Y", "Y")))

# F12: Linear regression models of voters’ preferences for cuts to government spending, dropping one's own party as control
m1.0 <- lm(cut_spend ~   l_ip_per414_mean * l_ip_per414_consensus + country + year + loggdp + unemp + growth + age + gender, weights = weight, data = df)
m1 <-   lm(cut_spend ~   l_ip_per414_mean * l_ip_per414_consensus + country + year + loggdp + unemp + growth + age + gender, weights = weight, data = subset(df, inc.group==0))
m2 <-   lm(cut_spend ~   l_ip_per414_mean * l_ip_per414_consensus + country + year + loggdp + unemp + growth + age + gender, weights = weight, data = subset(df, inc.group==1))
m3 <-   lm(cut_spend ~   l_ip_per414_mean * l_ip_per414_consensus + country + year + loggdp + unemp + growth + age + gender, weights = weight, data = subset(df, inc.group==2))

## Robust SE: Baseline model
r1.0 <-  sqrt(diag(vcovHC(m1.0, type = "HC1")))
r1 <-  sqrt(diag(vcovHC(m1, type = "HC1")))
r2 <-  sqrt(diag(vcovHC(m2, type = "HC1")))
r3 <-  sqrt(diag(vcovHC(m3, type = "HC1")))

# Regression
stargazer(m1.0, m1, m2, m3,
          se = list(r1.0, r1, r2, r3),
          keep.stat=c("n", "adj.rsq"),
          type = "text", 
          # out = "../05_figures/v3/appendix_f12.html",
          omit = c("country", "year"),
          covariate.labels = c("Economic orthodoxy (average, t-1)", "Economic orthodoxy (consensus, t-1)", "GDP (logged) (t)", "Unemployment rate (t)", "Growth (t)", "Age", "Female", "Economic orthodoxy average * consensus"),
          dep.var.labels = c("Preference for spending cuts"),
          add.lines=list(c("Income group", "All",  "L", "M", "H"),
                         c("Country-years", "60", "60", "60", "60"),
                         c("Fixed effects", "C&Y", "C&Y", "C&Y", "C&Y", "C&Y", "C&Y", "C&Y", "C&Y", "C&Y")))

# F13: Linear regression models of voters’ preferences for spending in pensions, health and unemployment, dropping one's own party as control
m1.0 <- lm(spend_polarizing ~   l_ip_per414_mean * l_ip_per414_consensus + country + year + loggdp + unemp + growth + age + gender, weights = weight, data = df)
m1 <-   lm(spend_polarizing ~   l_ip_per414_mean * l_ip_per414_consensus + country + year + loggdp + unemp + growth + age + gender, weights = weight, data = subset(df, inc.group==0))
m2 <-   lm(spend_polarizing ~   l_ip_per414_mean * l_ip_per414_consensus + country + year + loggdp + unemp + growth + age + gender, weights = weight, data = subset(df, inc.group==1))
m3 <-   lm(spend_polarizing ~   l_ip_per414_mean * l_ip_per414_consensus + country + year + loggdp + unemp + growth + age + gender, weights = weight, data = subset(df, inc.group==2))

## Robust SE: Baseline model
r1.0 <-  sqrt(diag(vcovHC(m1.0, type = "HC1")))
r1 <-  sqrt(diag(vcovHC(m1, type = "HC1")))
r2 <-  sqrt(diag(vcovHC(m2, type = "HC1")))
r3 <-  sqrt(diag(vcovHC(m3, type = "HC1")))

# Regression
stargazer(m1.0, m1, m2, m3,
          se = list(r1.0, r1, r2, r3),
          keep.stat=c("n", "adj.rsq"),
          type = "text", 
          # out = "../05_figures/v3/appendix_f13.html",
          omit = c("country", "year"),
          covariate.labels = c("Economic orthodoxy (average, t-1)", "Economic orthodoxy (consensus, t-1)", "GDP (logged) (t)", "Unemployment rate (t)", "Growth (t)", "Age", "Female", "Economic orthodoxy average * consensus"),
          dep.var.labels = c("Preference for spending in pensions, health and unemployment"),
          add.lines=list(c("Income group", "All",  "L", "M", "H"),
                         c("Country-years", "60", "60", "60", "60"),
                         c("Fixed effects", "C&Y", "C&Y", "C&Y", "C&Y", "C&Y", "C&Y", "C&Y", "C&Y", "C&Y")))